D-place FARM documentation: Module 3

library(png)

format(Sys.time(), format="%d_%b_%Y")

[1] "27_Feb_2017"

## First consolidate the available files into a single table
    
      path <- "~/Box Sync/Four model compare/Module 2 extinct"
           
     
           setwd(path)

Error in setwd(path) : cannot change working directory

setwd("~/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries")

The working directory was changed to /Users/Ty/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.

details <- file.info(list.files())
trimmed_details <- details[which(list.files() == list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]

[1] "Four_model_compare_results_extinct_27_Feb_2017.Rdata"

load(rownames(trimmed_details[ord,])[1])
extinct <- Concatenated_data
trimmed_details <- details[which(list.files() != list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]

[1] "Four_model_compare_results_27_Feb_2017_crop_to_2861.Rdata"

load(rownames(trimmed_details[ord,])[1])
extant <- Concatenated_data

for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
    extinct[which(is.nan(as.numeric(as.character(extinct[, i]))) == TRUE), i] <- NA
}
for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
    extant[which(is.nan(as.numeric(as.character(extant[, i]))) == TRUE), i] <- NA
}
i <- 19
for(i in c(20,21,24,25,26,27)){
    extinct[which(as.numeric(as.character(extinct[, i])) == 0), i] <- NA
}
for(i in c(20,21,24,25,26,27)){
    extant[which(as.numeric(as.character(extant[, i])) == 0), i] <- NA
}
xlimit <- c(0,1)
ylimit <- c(0,600)
maincex <- 0.9
png(file="Global_success_rate_per_parameter.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(5,4), mar=c(3,3,3,0))
hist(as.numeric(as.character(extinct[,9])), main="speciation of F in F env", col=adjustcolor("firebrick", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,9])), main="speciation of F in F env", col=adjustcolor("cornflowerblue", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,10])), main="speciation of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[,10])), main="speciation of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,11])), main="speciation of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,11])), main="speciation of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,12])), main="speciation of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[,12])), main="speciation of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
#######
hist(as.numeric(as.character(extinct[, 14])), main="extinction of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 14])), main="extinction of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 15])), main="extinction of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 15])), main="extinction of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 16])), main="extinction of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 16])), main="extinction of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 17])), main="extinction of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 17])), main="extinction of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 29])), main="arisal of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 29])), main="arisal of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 30])), main="arisal of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 30])), main="arisal of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 31])), main="arisal of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 31])), main="arisal of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 32])), main="arisal of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 32])), main="arisal of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 20])), main="Diffusion: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 20])), main="Diffusion: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 21])), main="Diffusion: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 21])), main="Diffusion: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)

hist(as.numeric(as.character(extant[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
####
hist(as.numeric(as.character(extinct[, 24])), main="Takeover: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 24])), main="Takeover: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 25])), main="Takeover: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 25])), main="Takeover: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 26])), main="Takeover: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 26])), main="Takeover: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 27])), main="Takeover: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 27])), main="Takeover: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
dev.off()

null device 
          1 

png(file="extiction minus extant per outcome.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(3,1))
plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
plot(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))

plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
points(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6)
dev.off()

null device 
          1 

load('~/Box Sync/colliding ranges/Simulations_humans/Available trees/real.analysis.RData')
#Concatenated_data <- Concatenated_data[Concatenated_data[, 2] == "stats.no.bTO", ]
#Concatenated_data <- Concatenated_data[Concatenated_data[, 6] != "05", ]
# Concatenated_data[, 6] <- as.numeric(Concatenated_data[, 6])
# # Concatenated_data[original[, 2] == "background_takeover", 6] <-  Concatenated_data[original[, 2] == "background_takeover", 6] + 4
Concatenated_data[, 6] <- factor(Concatenated_data[, 6])
#head(Concatenated_data)
#names(Concatenated_data)
PCAdata <- Concatenated_data[, -(1:35)]
PCAdata <- PCAdata[, -12]
PCAdata <- apply(PCAdata, 2, as.numeric)
remove <- apply(is.na(PCAdata), 1, any)
PCAdata <- PCAdata[!remove, ]
# Predictions
library(randomForest)
data.analysis.comp2 <- data.frame("Model" = as.factor(Concatenated_data[!remove, 6]),
                                  PCAdata)
data.analysis.comp2$sprate <- data.analysis.comp2$trait_1_speciation/data.analysis.comp2$trait_2_speciation
data.analysis.comp2$extrate <- data.analysis.comp2$trait_1_extinction/data.analysis.comp2$trait_2_extinction
#load("Real_phy/real.analysis.RData")
a <- as.data.frame(real.analysis$results_summary_of_single_value_outputs)
a$sprate <- a$trait_1_speciation / a$trait_2_speciation
a$extrate <- a$trait_1_extinction / a$trait_2_extinction
data.analysis.comp3 <- data.analysis.comp2[, -c(2, 13:14, 16:20)]
#data.analysis.comp3 <- data.analysis.comp3[data.analysis.comp3$Model %in% 1:4, ]
#data.analysis.comp3$Model <- factor(data.analysis.comp3$Model)
#sub <- unlist(lapply(as.list(c(1:4)), function(x, y) {
#  sample(which(y$Model == x), min(table(data.analysis.comp3$Model)))},
#  y = data.analysis.comp3))
# data.analysis.comp3 <- data.analysis.comp3[sub, ]
fun <- function(x, y, per = .33) {sample(which(y$Model == x), round(table(y$Model)[1]*per))}
sub.test <- unlist(lapply(as.list(paste0(0, c(1:4))), fun,
                          y = data.analysis.comp3))
test2 <- data.analysis.comp3[sub.test, 2:ncol(data.analysis.comp3)]
test1 <- data.analysis.comp3[sub.test, 1]
train <- data.analysis.comp3[-sub.test, ]
(fit <- randomForest(Model ~ ., data=train, xtest = test2, ytest = test1, 
                    importance=TRUE, ntree=2400, keep.forest = TRUE, replace=TRUE, norm.votes=FALSE))

Call:
 randomForest(formula = Model ~ ., data = train, xtest = test2,      ytest = test1, importance = TRUE, ntree = 2400, keep.forest = TRUE,      replace = TRUE, norm.votes = FALSE) 
               Type of random forest: classification
                     Number of trees: 2400
No. of variables tried at each split: 4

        OOB estimate of  error rate: 36.25%
Confusion matrix:
     01   02   03 04 class.error
01 2796  579  473  8   0.2748963
02  669 2642  214 29   0.2566123
03  809  295 1060 25   0.5157606
04   63  441  121 55   0.9191176
                Test set error rate: 48.7%
Confusion matrix:
     01   02  03 04 class.error
01 1430  258 208  3   0.2469721
02  332 1440 105 22   0.2417062
03  679  258 939 23   0.5055292
04  183 1257 371 88   0.9536598

predictions <- predict(fit, 
                       a,
                       type="prob")

plot(fit, ylim=c(0,1))

labs <- c("Basic", "+Diffusion", "+Takeover", "+Diffusion +Takeover")

# bar plot
png("Prob_aus.png", width = 25, height = 25, res = 300, units = "in")
par(mar = c(8, 8, 1, 1))
pred <- setNames(as.numeric(predictions), labs)
cols <- rev(c("darkgreen", "red", "blue", "darkorange1"))
barplot(pred, col = cols, ylab = "Proability", cex.lab = 3, cex.names = 2)
dev.off()

null device 
          1 

prop

          01        02        03        04
01 62.280130 13.723696 21.829763  4.678007
02 13.289902 64.409881  8.627678 27.764277
03 20.065147  4.757548 53.503185 13.912515
04  4.364821 17.108875 16.039375 53.645200

importance(fit)

                                                   01        02        03
Pylo_diversity_is_sum_of_BL                  36.27731 26.313769 38.900228
average_phylogenetic_diversity_is_mean_of_BL 42.47717 32.482005 38.519904
variance_Pylo_diversity_is_variance_of_BL    42.15824 16.661635  9.266383
F_quadratic_entropy_is_sum_of_PD             24.95204 19.851036 36.048651
Mean_pairwise_distance                       31.03797 25.728248 33.157743
variance_pairwise_distance                   31.43706 29.867035 34.057818
Evolutionary_distinctiveness_sum             36.73711 26.605041 37.480353
mean_Phylogenetic_isolation                  40.52562 33.527197 41.801064
variance_Phylogenetic_isolation              24.57867  8.325148 37.260442
gamma                                        26.36037 21.591877 63.126688
extinction_per_speciation                     7.99549 41.950330 54.475002
transition_from_trait_1_to_2                 22.53100 22.856490 41.900317
transition_from_trait_2_to_1                 33.37167 33.715083 45.090491
transition_rate_ratio_1to2_over_2to1         50.20133 36.600927 35.430284
Phylogenetic_signal                          64.80645 43.118360 53.232738
spatial.tests.fora                           43.20987 86.600931 89.746230
spatial.tests.dom                            82.03166 30.318014 47.558336
prevalence                                   54.14112  9.616747 35.927323
sprate                                       41.34870 28.286586 85.062626
extrate                                      28.88581 24.078228 21.624649
                                                   04 MeanDecreaseAccuracy
Pylo_diversity_is_sum_of_BL                  45.76293             62.84990
average_phylogenetic_diversity_is_mean_of_BL 32.73694             67.29494
variance_Pylo_diversity_is_variance_of_BL    27.05678             54.35726
F_quadratic_entropy_is_sum_of_PD             43.18824             68.31908
Mean_pairwise_distance                       42.30806             63.22960
variance_pairwise_distance                   22.73712             60.82865
Evolutionary_distinctiveness_sum             44.98800             63.01409
mean_Phylogenetic_isolation                  33.38172             68.79250
variance_Phylogenetic_isolation              22.99055             49.94780
gamma                                        51.68400             89.28718
extinction_per_speciation                    54.09510             75.50858
transition_from_trait_1_to_2                 41.11630             58.91197
transition_from_trait_2_to_1                 40.76374             67.48847
transition_rate_ratio_1to2_over_2to1         42.40223             85.56456
Phylogenetic_signal                          40.17895             85.09066
spatial.tests.fora                           76.44656            125.15481
spatial.tests.dom                            53.97271            104.59956
prevalence                                   58.91447             79.90220
sprate                                       36.41170             92.93119
extrate                                      36.77120             58.20976
                                             MeanDecreaseGini
Pylo_diversity_is_sum_of_BL                          210.2330
average_phylogenetic_diversity_is_mean_of_BL         182.1725
variance_Pylo_diversity_is_variance_of_BL            164.4793
F_quadratic_entropy_is_sum_of_PD                     185.9903
Mean_pairwise_distance                               206.3388
variance_pairwise_distance                           183.6735
Evolutionary_distinctiveness_sum                     207.9225
mean_Phylogenetic_isolation                          184.8331
variance_Phylogenetic_isolation                      167.5043
gamma                                                213.3865
extinction_per_speciation                            251.9260
transition_from_trait_1_to_2                         255.1930
transition_from_trait_2_to_1                         262.6217
transition_rate_ratio_1to2_over_2to1                 211.2503
Phylogenetic_signal                                  303.1266
spatial.tests.fora                                   395.4726
spatial.tests.dom                                    250.7786
prevalence                                           226.0943
sprate                                               244.7882
extrate                                              190.3293

# Variables importance
imp <- importance(fit)
imp <- apply(imp, 2, function(x) (x - min(x))/(max(x) - min(x)))
imp <- imp[sort(imp[, 5], index.return = TRUE, decreasing = TRUE)$ix, ]
names <- rownames(imp)
names[names == "spatial.tests.fora"] <- "Space F"
names[names == "spatial.tests.dom"] <- "Space D"
names[names == "sprate"] <- "Sp(ratio)"
names[names == "transition_from_trait_1_to_2"] <- "TR(1-2)"
names[names == "transition_from_trait_2_to_1"] <- "TR(2-1)"
names[names == "Phylogenetic_signal"] <- "PhySig(D)"
names[names == "Evolutionary_distinctiveness_sum"] <- "EDsum"
names[names == "Pylo_diversity_is_sum_of_BL"] <- "PDsum"
names[names == "transition_rate_ratio_1to2_over_2to1"] <- "TR(ratio)"
names[names == "gamma"] <- "Gamma"
names[names == "mean_Phylogenetic_isolation"] <- "MPI"
names[names == "extrate"] <- "Ext(ratio)"
names[names == "average_phylogenetic_diversity_is_mean_of_BL"] <- "PDmean"
names[names == "extinction_per_speciation"] <- "DR"
names[names == "variance_Phylogenetic_isolation"] <- "VPI"
names[names == "F_quadratic_entropy_is_sum_of_PD"] <- "F"
names[names == "Mean_pairwise_distance"] <- "MPD"
names[names == "variance_Pylo_diversity_is_variance_of_BL"] <- "PDvar"
names[names == "variance_pairwise_distance"] <- "VPD"
png("var_import_all.png", width = 25, height = 25, unit="in", res=300)
par(mar = c(10, 18, 1, 1))
plot(x = rev(imp[, 5]), y = 1:nrow(imp), type = "l", yaxt = "n", 
     ylab = "", xlab = "Variable Importance",
     xlim = c(0, 1), lwd = 2, cex.lab = 4)
for (i in 1:nrow(imp)) {
  abline(h = i, lty = 3, col = "gray80")
}
abline(v = seq(0, 1, 1/19), lty = 3, col = "gray80")
lines(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", lwd = 2)
lines(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", lwd = 2)
lines(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", lwd = 2)
lines(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", lwd = 2)
lines(x = rev(imp[, 5]), y = 1:nrow(imp), lwd = 3)
points(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", cex = 2)
points(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", cex = 2)
points(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", cex = 2)
points(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", cex = 2)
points(x = rev(imp[, 5]), y = 1:nrow(imp), pch = 20, cex = 3)
text(y = 1:nrow(imp), x = par("usr")[1] - .17, labels = rev(names),
     srt = 0, pos = 4, xpd = T, cex = 4)
dev.off()

null device 
          1 

par(mfrow=c(2,3))
# Box plots
boxplot(spatial.tests.fora ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(spatial.tests.dom ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(log(sprate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(log(extrate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$extrate), col = "red", lty = 2)

boxplot(log(transition_rate_ratio_1to2_over_2to1) ~ Model, data = data.analysis.comp3)
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(Phylogenetic_signal ~ Model, data = data.analysis.comp3, ylim = c(0, 1))
abline(h = a$Phylogenetic_signal, col = "red", lty = 2)

str(fit)

List of 19
 $ call           : language randomForest(formula = Model ~ ., data = train, xtest = test2, ytest = test1,      importance = TRUE, ntree = 2400, keep.forest = TRUE)
 $ type           : chr "classification"
 $ predicted      : Factor w/ 4 levels "01","02","03",..: 1 3 1 1 3 3 3 1 3 1 ...
  ..- attr(*, "names")= chr [1:8885] "1" "2" "3" "4" ...
 $ err.rate       : num [1:2400, 1:5] 0.547 0.547 0.549 0.535 0.533 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:5] "OOB" "01" "02" "03" ...
 $ confusion      : num [1:4, 1:5] 1392 298 434 90 230 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. ..$ : chr [1:5] "01" "02" "03" "04" ...
 $ votes          : matrix [1:8885, 1:4] 0.549 0.346 0.906 0.902 0.22 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:8885] "1" "2" "3" "4" ...
  .. ..$ : chr [1:4] "01" "02" "03" "04"
  ..- attr(*, "class")= chr [1:2] "matrix" "votes"
 $ oob.times      : num [1:8885] 907 876 855 911 850 863 879 836 883 902 ...
 $ classes        : chr [1:4] "01" "02" "03" "04"
 $ importance     : num [1:20, 1:6] 0.0248 0.0262 0.0149 0.012 0.0186 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..$ : chr [1:6] "01" "02" "03" "04" ...
 $ importanceSD   : num [1:20, 1:5] 0.000668 0.000584 0.00033 0.000355 0.000494 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..$ : chr [1:5] "01" "02" "03" "04" ...
 $ localImportance: NULL
 $ proximity      : NULL
 $ ntree          : num 2400
 $ mtry           : num 4
 $ forest         :List of 14
  ..$ ndbigtree : int [1:2400] 3853 3761 3829 3805 3823 3977 3887 3917 3887 3907 ...
  ..$ nodestatus: int [1:4035, 1:2400] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ bestvar   : int [1:4035, 1:2400] 13 8 11 16 7 17 10 7 18 11 ...
  ..$ treemap   : int [1:4035, 1:2, 1:2400] 2 4 6 8 10 12 14 16 18 20 ...
  ..$ nodepred  : int [1:4035, 1:2400] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ xbestsplit: num [1:4035, 1:2400] 193.4137 0.0287 0.9995 25.8679 33.4238 ...
  ..$ pid       : num [1:4] 1 1 1 1
  ..$ cutoff    : num [1:4] 0.25 0.25 0.25 0.25
  ..$ ncat      : Named int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  ..$ maxcat    : int 1
  ..$ nrnodes   : int 4035
  ..$ ntree     : num 2400
  ..$ nclass    : int 4
  ..$ xlevels   :List of 20
  .. ..$ Pylo_diversity_is_sum_of_BL                 : num 0
  .. ..$ average_phylogenetic_diversity_is_mean_of_BL: num 0
  .. ..$ variance_Pylo_diversity_is_variance_of_BL   : num 0
  .. ..$ F_quadratic_entropy_is_sum_of_PD            : num 0
  .. ..$ Mean_pairwise_distance                      : num 0
  .. ..$ variance_pairwise_distance                  : num 0
  .. ..$ Evolutionary_distinctiveness_sum            : num 0
  .. ..$ mean_Phylogenetic_isolation                 : num 0
  .. ..$ variance_Phylogenetic_isolation             : num 0
  .. ..$ gamma                                       : num 0
  .. ..$ extinction_per_speciation                   : num 0
  .. ..$ transition_from_trait_1_to_2                : num 0
  .. ..$ transition_from_trait_2_to_1                : num 0
  .. ..$ transition_rate_ratio_1to2_over_2to1        : num 0
  .. ..$ Phylogenetic_signal                         : num 0
  .. ..$ spatial.tests.fora                          : num 0
  .. ..$ spatial.tests.dom                           : num 0
  .. ..$ prevalence                                  : num 0
  .. ..$ sprate                                      : num 0
  .. ..$ extrate                                     : num 0
 $ y              : Factor w/ 4 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "names")= chr [1:8885] "1" "2" "3" "4" ...
 $ test           :List of 5
  ..$ predicted: Factor w/ 4 levels "01","02","03",..: 1 3 1 1 1 3 1 3 2 1 ...
  .. ..- attr(*, "names")= chr [1:188] "1410" "86" "2307" "1834" ...
  ..$ err.rate : num [1:2400, 1:5] 0.559 0.537 0.473 0.479 0.473 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:5] "Test" "01" "02" "03" ...
  ..$ confusion: num [1:4, 1:5] 27 6 10 4 6 24 1 8 11 2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. .. ..$ : chr [1:5] "01" "02" "03" "04" ...
  ..$ votes    : matrix [1:188, 1:4] 0.94 0.354 0.532 0.543 0.525 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:188] "1410" "86" "2307" "1834" ...
  .. .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. ..- attr(*, "class")= chr [1:2] "matrix" "votes"
  ..$ proximity: NULL
 $ inbag          : NULL
 $ terms          :Classes 'terms', 'formula'  language Model ~ Pylo_diversity_is_sum_of_BL + average_phylogenetic_diversity_is_mean_of_BL +      variance_Pylo_diversity_is_variance_of_BL + F_quadratic_entropy_is_sum_of_PD +  ...
  .. ..- attr(*, "variables")= language list(Model, Pylo_diversity_is_sum_of_BL, average_phylogenetic_diversity_is_mean_of_BL,      variance_Pylo_diversity_is_variance_of_BL, F_quadratic_entropy_is_sum_of_PD,  ...
  .. ..- attr(*, "factors")= int [1:21, 1:20] 0 1 0 0 0 0 0 0 0 0 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:21] "Model" "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" ...
  .. .. .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..- attr(*, "term.labels")= chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..- attr(*, "order")= int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "intercept")= num 0
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(Model, Pylo_diversity_is_sum_of_BL, average_phylogenetic_diversity_is_mean_of_BL,      variance_Pylo_diversity_is_variance_of_BL, F_quadratic_entropy_is_sum_of_PD,  ...
  .. ..- attr(*, "dataClasses")= Named chr [1:21] "factor" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:21] "Model" "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" ...
 - attr(*, "class")= chr [1:2] "randomForest.formula" "randomForest"